
library(latex2exp)
library(repr)
library(digest)
library(gridExtra)
library(mltools)
library(cowplot)
library(infer)
library(AER)
install.packages("glmtoolbox")
library(glmtoolbox)
suppressWarnings(suppressMessages(library(tidyverse)))
suppressWarnings(suppressMessages(library(dplyr)))
suppressWarnings(suppressMessages(library(broom)))
suppressWarnings(suppressMessages(library(GGally)))
suppressWarnings(suppressMessages(library(leaps)))
suppressWarnings(suppressMessages(library(glmnet)))
suppressWarnings(suppressMessages(library(faraway)))
suppressWarnings(suppressMessages(library(base)))

Loading required package: car

Loading required package: carData

Loading required package: lmtest

Loading required package: zoo


Attaching package: ‘zoo’


The following objects are masked from ‘package:base’:

    as.Date, as.Date.numeric


Loading required package: sandwich

Loading required package: survival

Updating HTML index of packages in '.Library'

Making 'packages.html' ...
 done


insurance <- read_csv("https://raw.githubusercontent.com/LukeXboy/STAT-301-Group-Project/main/insurance.csv")
head(insurance)

Rows: 1338 Columns: 7
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (3): sex, smoker, region
dbl (4): age, bmi, children, charges

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.


#Produces further information about the feature values within the dataset
summary(insurance)

      age            sex                 bmi           children    
 Min.   :18.00   Length:1338        Min.   :15.96   Min.   :0.000  
 1st Qu.:27.00   Class :character   1st Qu.:26.30   1st Qu.:0.000  
 Median :39.00   Mode  :character   Median :30.40   Median :1.000  
 Mean   :39.21                      Mean   :30.66   Mean   :1.095  
 3rd Qu.:51.00                      3rd Qu.:34.69   3rd Qu.:2.000  
 Max.   :64.00                      Max.   :53.13   Max.   :5.000  
    smoker             region             charges     
 Length:1338        Length:1338        Min.   : 1122  
 Class :character   Class :character   1st Qu.: 4740  
 Mode  :character   Mode  :character   Median : 9382  
                                       Mean   :13270  
                                       3rd Qu.:16640  
                                       Max.   :63770


#Combining similar levels of region for future feature selection
insurance_factor <- insurance %>%
    mutate(region = ifelse(insurance$region == "northwest" | insurance$region == "northeast", "north", "south"))
head(insurance_factor)

#Mutating the features to be categorical variables for further analysis
insurance_factor <- insurance_factor %>%
    mutate(sex = as.factor(sex)) %>%
    mutate(smoker = as.factor(smoker)) %>%
    mutate(region = as.factor(region)) 

#To check the feature values after mutating the variables
summary(insurance_factor)

      age            sex           bmi           children     smoker    
 Min.   :18.00   female:662   Min.   :15.96   Min.   :0.000   no :1064  
 1st Qu.:27.00   male  :676   1st Qu.:26.30   1st Qu.:0.000   yes: 274  
 Median :39.00                Median :30.40   Median :1.000             
 Mean   :39.21                Mean   :30.66   Mean   :1.095             
 3rd Qu.:51.00                3rd Qu.:34.69   3rd Qu.:2.000             
 Max.   :64.00                Max.   :53.13   Max.   :5.000             
   region       charges     
 north:649   Min.   : 1122  
 south:689   1st Qu.: 4740  
             Median : 9382  
             Mean   :13270  
             3rd Qu.:16640  
             Max.   :63770


#Implementing changes
insurance <- insurance_factor

#Printing number of rows related to the database
cat("Number of instances in the insurance dataset: ", nrow(insurance), "\n")
cat("Number of instances with missing values in the dataset: ", sum(is.na(insurance)), "\n")

Number of instances in the insurance dataset:  1338 
Number of instances with missing values in the dataset:  0


#Historgram to visualize response variable's, charges, distribution
options(repr.plot.width = 10, repr.plot.height = 8)
charges_hist <- insurance %>%
    ggplot(aes(x = charges)) + 
    geom_histogram(bins = 25, boundary = 0, colour = 'white') +
    labs(title = "Training Set Distribution of Patient's Medical Costs", 
         x = "Charges (USD)") +
    theme(text = element_text(size = 16)) +
    geom_vline(xintercept = mean(insurance$charges), colour = 'red', size = 3) 
charges_hist


#Comparing correlation across all variables with a GGPairs plot
options(repr.plot.width = 15, repr.plot.height = 15)
all_corr_plot <- insurance %>%
    ggpairs() +
  labs(title = "Distribution and Correlation across Variables",
       x = "Variable",
       y = "Variable") +
  theme(
    text = element_text(size = 15),
    plot.title = element_text(face = "bold"),
    axis.title = element_text(face = "bold"),
    axis.text.x = element_text(angle = 70, vjust = 1, hjust=1) #To ensure labels do not overlap
  )
all_corr_plot

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.


#Splitting the data by labelling the data by ID
set.seed(123456)
insurance$ID <- rownames(insurance)
insurance_training <- sample_n(insurance, size = nrow(insurance) * 0.70, replace = FALSE)
insurance_testing <- anti_join(insurance, insurance_training, by = "ID")

insurance_training <- select(insurance_training, -ID)
insurance_testing <- select(insurance_testing, -ID)

head(insurance_training)
head(insurance_testing)


#Printing number of rows within split dataset
cat("Number of instances in the training set: ", nrow(insurance_training), "\n")
cat("Number of instances in the testing set: ", nrow(insurance_testing), "\n")

Number of instances in the training set:  936 
Number of instances in the testing set:  402


#Estimate additive MLR model using all variables in the dataset
insurance_full_MLR <- lm(charges ~ ., data = insurance_training)
tidy(insurance_full_MLR)


#Obtaining out-of-sample predictions with the full model
insurance_test_pred_full_MLR <- predict(insurance_full_MLR, newdata = insurance_testing)
head(insurance_test_pred_full_MLR)


#Storing model evaluations
insurance_RMSE_models <- tibble(
    Model = "Full Additive Regression",
    R_MSE = rmse(
        preds = insurance_test_pred_full_MLR,
        actuals = insurance_testing$charges)
    )
insurance_RMSE_models


#Check if there is multicollinearity within the dataset using General Variance Inflation Factors
vif_insurance <- gvif(insurance_full_MLR)
round(vif_insurance, 3)

           GVIF df GVIF^(1/(2*df))
age      1.0143  1          1.0071
sex      1.0062  1          1.0031
bmi      1.0625  1          1.0308
children 1.0018  1          1.0009
smoker   1.0053  1          1.0026
region   1.0511  1          1.0252


#Using the forward selection algorithm to find a reduced LR model
insurance_forward_sel <- regsubsets(
    x = charges ~ ., nvmax = 6, 
    data = insurance_training,
    method = "forward",
    )

insurance_forward_summary <- summary(insurance_forward_sel)
insurance_forward_summary

Subset selection object
Call: regsubsets.formula(x = charges ~ ., nvmax = 6, data = insurance_training, 
    method = "forward", )
6 Variables  (and intercept)
            Forced in Forced out
age             FALSE      FALSE
sexmale         FALSE      FALSE
bmi             FALSE      FALSE
children        FALSE      FALSE
smokeryes       FALSE      FALSE
regionsouth     FALSE      FALSE
1 subsets of each size up to 6
Selection Algorithm: forward
         age sexmale bmi children smokeryes regionsouth
1  ( 1 ) " " " "     " " " "      "*"       " "        
2  ( 1 ) "*" " "     " " " "      "*"       " "        
3  ( 1 ) "*" " "     "*" " "      "*"       " "        
4  ( 1 ) "*" " "     "*" "*"      "*"       " "        
5  ( 1 ) "*" " "     "*" "*"      "*"       "*"        
6  ( 1 ) "*" "*"     "*" "*"      "*"       "*"


#Summary table of measures for all forward selection models
insurance_forward_summary_df <- tibble(
    n_input_variables = 1:6,
    RSS = insurance_forward_summary$rss,
    ADJ.R2 = insurance_forward_summary$adjr2,
    BIC = insurance_forward_summary$bic,
    Cp = insurance_forward_summary$cp
    )
insurance_forward_summary_df


#Illustrates Cp of forward selection models
options(repr.plot.width = 8, repr.plot.height = 8)
plot(summary(insurance_forward_sel)$cp,
  main = "Cp for forward selection",
  xlab = "Number of Input Variables", ylab = "Cp", type = "b", pch = 19,
  col = "red"
)


#Estimation with selected model from forward selection
insurance_reduced_MLR <- lm(charges ~ age + bmi + children + smoker + region, data = insurance_training
)

#Out-of-sample predictions using the reduced model 
insurance_test_pred_reduced_MLR <- predict(insurance_reduced_MLR, newdata = insurance_testing)
head(insurance_test_pred_reduced_MLR)


#Combining reduced model's results
insurance_RMSE_models <- rbind(
    insurance_RMSE_models,
    tibble(
    Model = "Reduced Forward Selection Additive Regression",
    R_MSE = rmse(
        preds = insurance_test_pred_reduced_MLR,
        actuals = insurance_testing$charges)
    )
    )
insurance_RMSE_models


#Obtaining coefficents of the reduced model
tidy(summary(insurance_reduced_MLR))

term	estimate	std.error	statistic	p.value
<chr>	<dbl>	<dbl>	<dbl>	<dbl>
(Intercept)	-10653.6330	1151.96202	-9.2482502	1.524723e-19
age	241.0604	14.44958	16.6828650	7.448097e-55
sexmale	-313.8132	402.86313	-0.7789573	4.362032e-01
bmi	317.5317	33.42351	9.5002491	1.719942e-20
children	513.3964	168.26007	3.0512073	2.343990e-03
smokeryes	23962.2000	504.52291	47.4947711	9.059586e-251
regionsouth	-1004.6347	412.39565	-2.4360943	1.503383e-02

n_input_variables	RSS	ADJ.R2	BIC	Cp
<int>	<dbl>	<dbl>	<dbl>	<dbl>
1	50967889939	0.6233669	-901.3072	418.680882
2	38792983690	0.7130275	-1149.9528	98.038270
3	35659260671	0.7359263	-1221.9509	16.992655
4	35301815892	0.7382926	-1224.5390	9.520146
5	35078675062	0.7396672	-1223.6326	5.606774
6	35055778451	0.7395570	-1217.4021	7.000000

term	estimate	std.error	statistic	p.value
<chr>	<dbl>	<dbl>	<dbl>	<dbl>
(Intercept)	-10793.5098	1137.63963	-9.487635	1.916535e-20
age	241.3615	14.44135	16.713221	4.956678e-55
bmi	316.6875	33.39887	9.481981	2.013649e-20
children	512.6034	168.22142	3.047195	2.375121e-03
smokeryes	23935.7995	503.27681	47.559909	2.824810e-251
regionsouth	-1002.8241	412.30190	-2.432257	1.519283e-02

The Personal Attributes of A Patient that Influences One's Medical Costs¶

Introduction¶

Preliminary Results¶

Importing Modules¶

Downloading Dataset¶

Feature Descriptions and Wrangling the Dataset¶

Variable Distributions and Correlations¶

Splitting Data into Training and Testing Sets¶

Methods: Plan¶

Results¶

Estimating Multilinear Regression Model¶

Mutlicollinearity¶

A Smaller Model using Forward Selection¶

Methods: Results¶

Discussion¶

References¶

age	sex	bmi	children	smoker	region	charges
<dbl>	<fct>	<dbl>	<dbl>	<fct>	<fct>	<dbl>
32	male	31.50	1	no	south	4076.497
59	male	27.50	1	no	south	12333.828
42	female	25.30	1	no	south	7045.499
42	male	34.10	0	no	south	5979.731
18	male	35.20	1	no	south	1727.540
18	female	26.73	0	no	south	1615.767

	GVIF	df	GVIF^(1/(2*df))
age	1.014	1	1.007
sex	1.006	1	1.003
bmi	1.062	1	1.031
children	1.002	1	1.001
smoker	1.005	1	1.003
region	1.051	1	1.025

Model	R_MSE
<chr>	<dbl>
Full Additive Regression	5888.656
Reduced Forward Selection Additive Regression	5882.354